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iD ' Abstract 






The problem of finding the optimal tapering of a free (non-supported) javelin is 
described and solved. For the optimal javelin, the lowest mode of vibration has the 
highest possible frequency. With this tapering inner damping will lead to the cessation 
r^ ' of the vibration at the fastest possible rate. The javelin is modeled as a beam of 

C^ , uniform material. The differential equations governing the vibration and the tapering 

of the beam are derived. These equations have a difficult singularity at the tips of the 
beam. A procedure using a similarity solution, as in 0], is used to solve this singular 
system, and the solution is found. The maximal frequency is found to be almost 5 times 
►^ , larger than the frequency of a cylindrical rod. 

^^ ' Keywords: vibrating beam, eigenvalue optimization, singular ODE, similarity solution, sta- 

^-^ . ble manifold. 

t^ ! 1 Introduction 

O 

The interest in the optimal design of columns, beams and plates has existed for many years. 
^ ' Euler started the rigorous study of the buckling load of columns and introduced the problem 

of designing the strongest column. Keller, in 1960 solved this problem [5] . In 1964 Keller 
and Niordson found the design of the tallest self- weighted column [^ . Others have continued 
studying the various qualities of bending rods and plates under various conditions. 

In this paper we find the optimal design of a non-supported beam (picture an Olympic 
javelin in mid-air). The aim is to find the design whose lowest mode of vibration has the 
largest frequency. The optimal design is shown to have a frequency that is greater than 
that of a constant cross-section beam by a factor of 5. 

To simplify the problem we make several working assumptions on the permissible designs of 
the column. The cross-sectional shapes at different points along the beam are assumed to 
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be geometrically similar with fixed orientation (see figure [T|). Furthermore, we assume that 
the cross-sectional shape is convex. The cross sectional area is allowed to vary throughout 
the length of the beam ("tapering"). While maximizing the frequency, we hold the total 
volume of the beam fixed. 

Working within linearized theory, it is sufficient to consider standing waves confined to a 
single plane. These standing waves and their temporal frequencies are solutions of an ODE 
eigenvalue problem. The frequencies are functionals of the beam shape. This analysis seeks 
the tapering of a beam with fixed length and volume, which maximizes the lowest frequency. 
Formally this is done by requiring that the frequency be stationary with respect to variation 
of the beam tapering. This gives an additional ODE which relates the tapering and the 
standing wave amplitude. 

The ODE's and boundary conditions form a closed system for the tapering, standing wave 
amplitude and frequency of the optimal beam. Unfortunately, they are difficult to solve. 
Naive shooting methods fail to get close to the end of the beam and therefore do not 
allow for corrections of the initial conditions to be made. More sophisticated boundary 
value problem solvers also fail to converge. In [8] Niordson solved a similar problem by 
converting the ODE to integral form and then performing an iteration which converges to 
the solution. This paper follows Niordson's paper loosely but since the boundary conditions 
(BC) are different and the method of solution is different, we present the full derivation and 
solution here. Having different BC means that although this problem has the same ODE's, 
the singularities at the tips are more severe in this case. We use the same method of solution 
shown before in [3]. First, in Section [2] the equations that characterize the optimal beam 
and the shape of vibration mode are found. As mentioned, these equations are nonlinear 
and singular at the tips of the beam. In Section [3] we reduce these singular equations to 
a regular system of ODE's that can be easily solved using standard numerical methods. 
A similarity solution to the equations is found and used to "peel away" the singularity 
at the tips. The resulting ODE's have a critical point and by starting near the critical 
point on its stable manifold, the equations are solved backwards numerically until the BC 
are satisfied. Since the stable manifold is two dimensional, a simple 1-parameter shooting 
algorithm employing a standard ODE integrator will determine the solution. 



2 Derivation of the Boundary Value Problem 

2.1 Setup 

Consider all possible beams, all of the same length and volume which are suspended without 
gravity or other external forces. The beams have various modes of vibration. What is the 
design of the beam whose first vibration mode has the largest frequency? To simplify, we 
solve the problem only for a specific class of permissible designs. We assume that the beam 
is thin (i.e. the characteristic width is much less than the length of the beam) and made 



of a homogeneous material. In addition, we only permit beams with geometrically similar, 
equally oriented and convex cross-sections. Lastly, we concern ourselves with the bending 
of the beam in a specified plane only. 

The beam is parameterized by arclength s, measured from one of the tips along the beam's 
center axis. The design information is contained in a single function a{s), the cross-sectional 
area of the beam at point s (see Figure [1]). At rest and under no stress, the center of beam 




Figure 1: The beam is assumed to be made of a homogeneous material with convex cross- 
section. The area of the cross-section at point s is denoted by a{s). The cross-sections are 
geometrically similar and equally oriented. 

is assumed to lie on the x— axis. The beam configuration at time t is specified by u{s,t), 
measuring the vertical displacement of the point s from the x— axis (see Figure [2]). 

First, we think of the cross sectional area, a{s), as given. The total volume of the beam is 
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where L is the total length of the beam. 



a{s) ds, 



(1) 



2.2 Lagrangian 

A beam design given by a{s) determines the vibration modes of the beam. To find the ODE 
that governs the vibration, we write the Lagrangian, given by the difference between the 
kinetic and potential energy of the beam: 
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Hs)u'is + 7;P(^is)uf \ ds- 



(2) 




Figure 2: A "snapshot" of the vibrating beam at time t. The vertical displacement of the 
beam is given by u{s,t). The position along the beam is parameterized by arclength s, 
measured from one of the tips of the beam. 



Here, p is the mass density of the material. The function b{s) is the bending modulus which 
is proportional to a^(s). Specifically, 



b{s) = cEa^{s), 



(3) 



where c is a dimensionless constant that depends on cross-section shape, and E is the 
Young's modulus of the material. Using separation of variables we write the deflection 
function u{s,t) as a product of a standing wave amplitude function, y{s) and cos{ujt), and 
average the Lagrangian ([2]) over a temporal period: 



^[y] 



W /" / I 1 --,2/„\„.2 „„„2/, ..\ , -'- „„/„\ „.2, ,2„- 2 



2tt 



Jo 



cEa {s)ygg cos (ujt) + - p a{s) y u sin (wt) ^ dsdt. (4) 



^j {-cEa'^is) yl + pa{s) y^uj^ } ds. 



(5) 



The average Lagrangian can be written using non-dimensional variables by implementing 
the units in the scaling table: 



Variable a s 



y 



jC 



Unit 



L 2 



2V ScEV^ 



The dimensionless versions of ([5]) and ([T|) are 

- 1 r^ 

^[y\ = lj {-ayL + A^ay^} ds, (6) 

V[a] = I ads = 2 (7) 
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Here, the square of the non-dimensional frequency is A^ = Cgg^y ■ 

2.3 Vibration ODE 

The Euler equations of ([6]) constitute a boundary value problem (BVP) for y{s): 

(a2y,,)^^-A2ay = 0, in < s < 2, (8) 

a^Vss = at s = 0, 2, {a^yss)s = at s = 0, 2. (9) 

Physically, the BC express the absence of torque and force at the ends. Although the 
Euler equations were derived from the average Lagrangian, finding the ODE from the full 
Lagrangian and then using separation of variables will lead to the same equations for y{s). 

Heuristically, it is reasonable to expect the optimal beam shape has a{s) even. We also 
expect the fundamental mode to be an even standing wave. This allows us to solve the 
problem on the interval < s < 1. The volume constraint ([7]) reduces to 

1 
a{s)ds = l. (10) 



At the endpoint s = 1 we impose symmetry boundary conditions on y(s), so the eigenvalue 
problem for the shape of standing waves is 

{a^yss)^^-\^ay = Q, 0<s<l (11) 

a^yss = 0, at s = 0, [a^yss)^ = 0, at s = 0, (12) 

y.(l) = 0, ysss{l) = 0. (13) 

2.4 The Frequency of a Cylindrical Javelin 

If the cross-section a{s) is a constant, the javelin is a simple cylinder. In this case we can 
solve the problem (almost) analytically. This will give us a reference frequency to compare 
with later. To find the frequency of a cylindrical javelin, it is more straightforward to shift 
the origin of s and solve on the interval [—1,1]. The even solutions to the ODE are 



y{s) = A COS (VXs) + B cosh (VXs) . (14) 

The BC yield a constraint on A: 

- tan (VX) = tanh (^^ (15) 

Solving this equation numerically for the smallest (nonzero) A gives 

A w 5.5933. (16) 

This is the non-dimensional frequency of the cylindrical javelin. The standing wave shape 
of the cylindrical javelin is shown in figure [3] 
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Figure 3: The shape of the first vibration mode of the cylindrical javelin. 



2.5 McLximizing the Frequency 

Up to this point, we found an eigenvalue problem that implicitly determines a functional, 
A [a]. For each cross-sectional area function a{s) it defines the frequency of a standing wave 
solution of a beam tapered according to it. Although the formula is implicit and we can by 
no means give an explicit formula for X[a], we would like to find the function a{s) so X[a] is 
stationary with respect to variations of a{s). To do this, we find the functional derivative 
of A [a] (with respect to a) and write 
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6a 



H. 



(17) 



In (J17p . V[a] is the volume functional (jlOp and ^u is a Lagrange multiplier associated with 
the volume constraint. To find 4^, we introduce a small variation to the cross-sectional area 
6a{s). Let 6\ and 6y be the resulting variations in A and y. Assuming the the resulting 
variations are small when 5a is small, the linear variational equations which follow from 
dTTHISD are 



{2a5ayss + OL 5yss)ss — '2X6Xay — X 6ay — X^a6y = 0, 



(18) 



2a6ayss+0' Si/ss = 0, at s = 0, (2a6ayss + a Si/ss) =0, at s = 0, (19) 

6ysil) = 0, 6y,ssil) = 0. (20) 

This is a linear, inhomogeneous BVP for 5y. Its solvability condition determines the re- 
lationship between 6a and SX and thus gives ^. The solvability condition is found by 
multiplying equation (jlSp by y and integrating from to 1. Integration by parts, use of the 
BC, and rearrangement give 

/ {{a'^yss)ss-X^ay} 6yds+ [ {2ayl - X^ y^] 5ads = 6X [ 2Xay^ds. (21) 
Jo Jo Jo 



The first integral in the left-hand integral vanishes due to (jlip . The remaining terms give 
the desired relationship between Sa and 6X: 
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(22) 



or 

5X ^ 2ayl-X^y^ 

5a 2Xf^ay^ds' 

Substituting the expression for ^ into equation (|17|) yields an integro-differential equation 
that characterizes the optimal tapering a(s), 

2ayl^-X'^y'^ = 2^iX ay^dr (24) 

Jo 

A short calculation shows that /x = A (see Appendix |A]). Since the RHS of (p^ is indepen- 
dent of s, this integral equation can be transformed into an ODE by differentiating it once 
with respect to s, 

2(ayL),-A' (y'), = 0. (25) 

It would seem that in order to remain equivalent to the integro-differential equation (j24p . 
an additional BC should be added. In fact, the volume constraint on a is enough. See 
Appendix [XI 

The volume constraint (jlOp . BVP (|llH13p and equation (|25p . characterize the tapering of 
the javelin with highest frequency. These equations are singular at the tip, s = 0, due to 
the BC (I12p and a direct numerical approach fails to give a solution. We will now find 
a similarity solution that will remove the singularity by transforming the ODE into an 
autonomous system which can be solved using a simple ODE solver. 



3 Solution of the BVP 

In order to manage the derivatives more easily, we introduce a new variable, 

Lp = o?yss- (26) 



This is the non-dimensional torque. Since the torque, 99, goes to zero at the tip (due to the 
BC), and we expect a{s) to go to zero as well, we look for an algebraic relation between 
the variables and the distance from the tip. The limit to be examined is s ^ 0. First we 
change the constraint J^ ads = 1 into a BVP. This can be achieved by adding a variable 

b{s): 

b{s) = [ a{s')ds'. (27) 

Jo 

In terms of the newly defined variables, the BVP is 

(^ - a^ yss = 0, (28) 

ipss-X^ay = 0, (29) 



(30) 
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6 = 0, 
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^s -- 


= 0, 


at s = 
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= 0, 


at s = 


= 1 


b = 


= 1, 


at s = 


= 1 



a = 0, (31) 

(32) 
(33) 
(34) 

In the equations above, BC (j32p is the translation of (J12p . The BC (I33p are, in essence, 
symmetry BC. They follow from (J13l25p . Lastly, (j34p are the BC needed for the volume 
constraint. They follow immediatly from (1101 I27p . 

3.1 Similarity Solution 



The ODE's ()28H3ip have a similarity solution which satisfies the BC at the tip of the beam. 
(For information on similarity solutions see, for example, [1] or [2].) To find it, we examine 
the scaling relations among the variables. Let A,B,P,Y and S be the "units" of a,b,(p,y 
and s respectively. A balance of "units" in the equations (|28H3ip gives the relationships 











P = A'YS-^ 










PS-^ = AY, 










P^A~^S~' = Y^S~\ 










BS~^ = A. 


This system 


has 


a two- 


-parameter 


' family of solutions: 

A = S\ 

P = YS*^. 



This leads us to look for a solution (j28H3ip of the form 

a{s) = ttQS , 
b{s) = 6os^ 

y{s) = yos^. 

In the equations above, the exponent p is unknown. 
Substituting these equations into (f28H3T]l yields 



^'o = ^, (35) 
5 

V^o =p(p- l)?/oao' (36) 

7 = p(p-l)(p + 6)(p + 5), (37) 

7 = 2/(p_l)2, (38) 



where 



\2 

7=-. (39) 

ao 

Equations ([37|) and ([38]) yield a polynomial equation for p, 

p{p-l){p + Q){p + f>) = 2p^{p-lf. (40) 

Since we are looking for a real frequency, A^ must be positive. Clearly, oq must also be 
positive since a{s) is an area. Therefore, 7 must also be positive. This rules out the two 
trivial solutions to (j40p . p = and p = 1. The two other solutions are —2, 15. 

The solution p = 15 gives a vanishing LHS for the integral equation (j24p (as s ^ 0). This is 
not possible unless the constant RHS is also zero. Since the RHS of ()24p is positive, p = 15 
is not a solution of interest. This leaves us with the single possible solution p = —2. This 
solution yields 7 = 72 and gives rise to the following similarity solution: 

a{s) = ^s\ (41) 

Ks) = 1^^ (42) 

^{s) = yo^s\ (43) 

y{s) = yos"^. (44) 

It is easy to check that (|41fl44p solves the ODE system for all s and satisfies the BC at 
the tip (.8 = 0). It is also easy to check that this solution does not satisfy the BC at the 
midpoint {s = 1). We now use this similarity solution to remove the singularity from the 
ODE's, simplifying the equations to a point where a numerical solution is possible. 
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3.2 Peeling away the Singularity 

To analyze the solution of the full BVP (|28H34p . we "peel away" the similarity solution. 
This is done by the transformation to the variables a, (3, $, C, defined by 

a{s) = a{s)a{s), (45) 

b{s) = b{s)P{s), (46) 

ip{s) = ip{s)^s), (47) 

yis) = yis)Cis). (48) 

Substituting these expressions into the BVP (p8HM|) results in a BVP for a,/3, $, and (. 
Since the resulting equations are homogeneous in s we use t = — In s as the independent 
variable. In this variable, the ODE for a{t), I3{t), $(t), and C(t) are the autonomous system 
(AS), 

6<^-a'^{D + 2){D + 3)C = 0, (49) 

{D - 4){D - 3)<^ - 12aC = 0, (50) 

(4 + D)^-2C(2 + D)C = 0, (51) 

(5 - D)/3 - 5a = 0. (52) 

Here D is the derivative with respect to the variable t. The BC for this system are: 

^(t)e-^* -^ 0, e-^*(Z) - 4)$(t) -^ 0, /?(i)e"^* ^0, as i ^ cx) (53) 

{D-A)a = 0, {D + 2)C = 0, -^ = T2"' at i = 0. (54) 

(a) (b) (c) 

What have we gained by all these manipulations? First, we notice that A is no longer part 
of the ODE. It only appears in the BC at t = 0. This greatly simplifies the solution of 
the BVP. Also, we notice that the singularity at s = has been removed. The boundary 
conditions do not cause the variables to vanish and there is no delicate balance of terms. 
The similarity solution of the original BVP ()41H44p is represented by the critical point 
(a,/9, <1>,C) = (1,1,1,1) = 1- Since the similarity solution satisfies the BC at the tip, we 
look for a solution that satisfies the BC at t=0 and converges to 1 as t ^ 00. This means 
that we are looking for a solution on the stable manifold of the fixed point 1. 

As is shown in the next section, the stable manifold is two-dimensional. On the other hand, 
the BC (I54h ) . ()54b ) define a surface of co-dimension 2. Thus, these surfaces are expected to 
have discrete points of intersection. These points, via the BC (|54b ) determine a particular 
value of A. In our case we will find exactly one point and hence one possible value for A. 
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3.3 Stable Manifold 



To find the tangent plane of the stable manifold of 1 , we linearize the ODE around 1 and 
search for solutions of the form 

y = yoe^*- 

The directions with "^{q) < are stable. The linearization of the AS (|49ll52p is, 



V 



-12 6 

-12 (Z?-4)(Z?-3) 

-3(4 + L>) 2(4 + 1?) 

-5 (5-D) 



-{2 + D){3 + D) \ 



-12 

-2(4 + D) 





/ 



/ 6a \ 

613 

\sc J 



0. 



(55) 



Here, {6a,6P,6^,6C) are deviations from 1. We look for solutions of this system in the 
form 

i6a,6p,6^,6C) = {6ao,6Po,S^o,SCo)e'''. 

Substitution into ([55]) yields 



/ -12 6 

-12 (g-4)(g-3) 

-3(4 + g) 2(4 + ^) 

-5 (5 - g) 



-i2 + q)i3 + q) \ f 6ao \ 



-12 

-2(4 + g) 





/ 



6(3o 
6-^0 



0. 



(56) 



This system has a nonzero solution for {6ao,6PQ, 6^0,6(0) when the matrix in (j56p is sin- 
gular. This happens for 6 values oi q : qi = 0, q2 = 1, qs = —4, 54 = 5, ^5 ss 6.3523, 
q^ ss —5.3523. The corresponding solutions for {6ao, 6f3o, 6^0, ^Co) are given in Tabled) 
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^2 


S3 


^4 


^5 


^6 


5ao 
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0.68568 


-11.019 


SPo 





5 


5 


1 


-2.5352 


-5.3220 


<5$o 


1 


4 


-27 





1 


1 


-JCo 


1 


-2 


135 





-0.028525 


17.529 



Table 1: The stable and unstable directions around the critical point of the AS (J49I I& 



The values for 5*5 and Sq are approximate. The only two stable solutions are {q3,Ss), 
(^6) Sq); therefore, the plane tangent to the stable manifold is spanned by the two vectors 
53 and Sq. The unstable direction ^2 is due to the similarity solution (see [3]), and the 
unstable direction §4 is due to (3 representing an integral constraint on the solution, not 
actually coupled to the ODE (See Appendix lA.li ) 

To find the numerical solution of the BVP, we start near the fixed point P, on the plane 
tangent to the stable manifold, and solve the AS ()49H52p backwards in t. The stopping 
condition is that both BC (I54bl and (I54bl are satisfied at the same t. Since the stable 
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manifold is two-dimensional, we have a one-parameter family of solutions each starting at a 
different direction on the manifold. We use the shooting method to find the initial direction, 
so the resulting solution for (a, (3, $, C) satisfies both BC at the same t. Since the system is 
autonomous, we redefine this t to be zero. Once stopped, the value of A will be determined 
from (|54b ) and then the full solution follows using 



4 Numerical Results 

Here is the actual mechanism of the shooting method: First, the direction in the stable 
manifold is defined using a parameter € [0, 2tt): 

v{9) = sm{e)S3 + cos{0)Sq. (57) 

Next, the AS is solved backwards in t starting from 

xo = (l,l,l,l)+ef(^), (58) 

where e is a small parameter determining how close to the fixed point to start the solution. 
A value of e = 0.001 was used in this numerical solution. As the ODE is 6th order, initial 
values for ^j and $t are needed. For this we used the derivatives of the similarity solution: 
Ct = eqdCo and $t = eq6^o. 

The AS is solved using Matlab ode solver ode45 using default tolerances. Plotting —At 
for which each of the two BC are satisfied (for each value of 9), gives Figure HI We see that 
for some values of 6 one or both of the BC are never satisfied, while for others a BC can 
be satisfied several times. The two BC are satisfied for the same —At for a single value of 
6, around 7r/6. Using the Matlab non-linear solver, f zero, the value of 9 where the two 
BC are satisfied at the same t is found, 6 ~ 5.753. The solver was given — vr/G as the initial 
guess for 6. For this 9 the two BC are satisfied at At = —2.0429. From the value of P at 
s = 1 the value for A is found: A ~ 27.073. We can compare this value of A to the value 
for the simple cylinder. The non-dimensional frequency for a cylindrical rod is 5.5933, and 
therefore the optimized rod vibrates almost 5 times faster than the cylindrical one. The 
tapering of the optimal javelin is shown in Figure [5] along with the shape of its standing 
wave. 

5 Discussion 

We have shown how to use the similarity solution to remove the singularity from the dif- 
ferential equations and find a solution that would otherwise require an iterative method. 
The variational equations were derived under the assumption that the spectrum of the dif- 
ferential operator ()llH13p is discrete and therefore the variation will have a meaning. Cox 
and McCarthy have shown (for example in [3] , [7j ) that this is not always the case and that 
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Figure 4: The —At at which each of the relevant BC are satisfied as a function of 9. The 
continuous hne is the value of —At for which a^ = 0, and the broken line is the value for 
which Us = 0. The plot shows that not every direction leads to a nicely behaved solution. In 
some directions the solution explodes before one or both of the BC are satisfied. In others 
the solution has two values of t for which a particular BC is satisfied. 

special treatment due to the existence of a continuous spectrum may be necessary. The 
existence of the continuous spectrum is due to the singularly tapered tips and therefore any 
minimal amount of rounding of the tips will eliminate the continuous spectrum. 

Another possible inaccuracy in the above derivation is due to the basic assumption that the 
defiection is small. The deflection y ends up having a singularity at the tips of the beam 
and therefore can only be small away from the tips. This means that the linearization is a 
crude estimate at the tips. In addition the curvature was taken to be equal to yss, this is 
only true when ys <^1. Again, this assumption breaks down near the tips where the slope, 
ys, tends to infinity. 

The extension of this analysis, to optimizing higher modes, is not obvious. The second 
mode is expected to be anti-symmetric and can be found using other BC at the middle of 
the beam (y = instead of y^ = 0). Higher modes may have singularities at internal points. 
To solve this "contact conditions" governing the internal singularities must be derived and 
used to connect between different parts of the solution. 
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(a) 



(b) 



Figure 5: The shape of the stiffest rod and the shape of the fundamental standing wave. 
Since the fundamental standing wave has an s~^ singularity at the tip, s^y is plotted instead 
of y. Both plots do not continue all the way to s = 0. This is because the solution was 
started a small distance away from the fixed point P. The solution can be easily continued 
near s = using the similarity solution. The dashed line in the figure on the left is the 
similarity solution a(s), the actual solution slowly leaves this solution as s increases. 

A Appendix 

Calculation of the Lagrange Multiplier /i 



To calculate the Lagrange multiplier in (j24p . we multiply the equation by a, integrate and 
use the volume constraint (1101) 



(59) 



2 / a^ ylgds — X^ I ay^ ds = ii\ j ay^ dr. 



/ ay'^ ds = ii\ I a 'I 
Jo Jo 

Two integration by parts (and use of the BC) yields 

2/ {a'^ yss) gg y ds - X^ ay'^ds = fi\ ay'^dr. 
Jo Jo Jo 

Using the ODE ([U]) we get 



2A^ 



ay ds — X / ay ds = fj,X ay dr. 



(60) 



(61) 



'0 Jo Jo 

Thus fi = X. By this calculation one can also "go back" from the differential equation ()25p 
to the integro-differential equation (j24p . Integrating (|25p once gives 

2/ a^yl,ds-X^I ay^ ds = C. (62) 

Jo Jo 

Here C is an unknown constant. Multiplying by a, integrating and using the volume con- 
straint (1101) recovers the constant C. 
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A.l An Unstable Direction 

Solution S4 in table ([1]) has a suspicious form. The eigenvalue is 5, the exponent of the 
similarity solution, and the eigenvector has components only in /3 direction. This is because 
the original ODE are invariant under a shift of b by an additive constant. Shifting 6 by e 
translates to the /3 variable: 

b 

360 . 

= 1 + .^.- 



360 

1 + e—iT-e 



5t 



So we see that there is an unstable direction about the critical point that makes (3 increase 
exponentially with constant 5. 
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